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Abstract The fungal infection called chytridiomycosis, 
caused by Batrachochytrium dendrobatidis (Bd), has given rise to 
dramatic declines or extinctions of many amphibian species 
around the world; however, in Asia, this disease has shown 
a low zoospore load or scant mortality. One potential reason 
for this may be that certain unique community structures 
of amphibian skin symbiosis contribute to the outcome 
of the disease; nevertheless, we know very little about the 
microbiota in this region. In this study, we used skin swabs 
of five sympatric amphibian species that have various 
habitat preferences in Lishui, Zhejiang Province, a place in 
southeastern China, to explore the skin bacterial communities 
by using 16S rRNA amplicon sequencing. We detected a 
total of 1020 OTUs, belonging to 17 phyla, among which 
Proteobacteria, Firmicutes, Actinobacteria and Bacteroidetes 
dominated all five host species. Enterobacteriaceae 
and Exiguobacteraceae and the genera Escherichia and 
Exiguobacterium belonging to these two families were 
identified as the most abundant taxa on our focal species. 
The alpha diversity was significantly lower on the terrestrial 
species, and also the highly enriched Proteobacteria was 
found on the terrestrial species, Rana zhenhaiensis, whereas 
Actinobacteria, Bacteroidetes, Cyanobacteria and Firmicutes 
were more abundant on aquatic species than on the terrestrial 
species. Our results suggest that both host species and habitat 
sites are important factors driving skin microbial diversity 
and composition and that amphibians in China may 
harbour unique skin bacterial communities. This study helps 
elucidate amphibian skin microbial ecology, and with further 
efforts, the specific mechanism of the interaction between Bd 


and host amphibians in China could be elucidated. 
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1. Introduction 


Vertebrates act as habitats for diverse and complex microbial 
communities that can contribute to host health or, conversely, 
be detrimental to the host. This holobiont (host plus symbionts) 
has aroused widespread concern for the past few years 
(Rosenberg and Zilber-Rosenberg, 2016), such as in studies on 
different parts of the human body, for example, the intestine 
(Brooks et al, 2018; Vitetta et al, 2018), mouth (Vogtmann et al, 
2019), and vagina (Bernabeu et al., 2019; Fettweis et al., 2019). 
In addition, the skin is the largest organ of the body, along 
with microbiome on it, is a complex and dynamic ecosystem 
and is given great importance in human studies (Chen et al., 
2018; Byrd et al., 2018). In recent years, the development of 
next generation sequencing technology allowing such studies 
about skin microbial symbionts has also been performed in 
animals, such as mammals (Hoffmann et al., 2014; Older et al., 
2017), birds (Roggenbuck et al., 2014; Engel et al., 2018), reptiles 
(Hyde et al, 2016; Allender et al., 2018; Walker et al, 2019), fish 
(Carda-Diéguez et al., 2017; Minniti et al, 2017) and, especially, 
amphibians (McKenzie et al., 2012; Kueneman et al, 2014; Jani et 
al, 2017; Medina et al, 2019). Amphibian skin is among the best- 
studied systems due to the notorious fungal infection known 
as chytridiomycosis caused by Batrachochytrium dendrobatidis 
(hereafter Bd) (Berger et al, 1998; Longcore et al., 1999; O'Hanlon 
et al, 2018). 

Batrachochytrium dendrobatidis has been linked to dramatic 
declines of more than 500 amphibian species around the 
world, amounting for 6.2% of the 8118 described amphibian 
species, 90 of which have been extinction (Scheele et al, 2019; 
AmphibiaWeb, 2020). An increasing number of studies have 
found that the skin microbiome of amphibians provides the 
first line of defence against Bd through the production of anti- 
fungal compounds such as violacein, indole-3-carboxaldehyde, 
2,4- diacetylphloroglucinol (DAPG) and prodigiosin (Brucker et 
al, 2008a, b; Harris et al., 2009; Woodhams et al., 2018; Madison 
et al., 2019). Previous studies on amphibian skin symbiosis 
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have mainly focused on Australia, Mesoamerica and America, 
where most amphibian decline events have occurred, while this 
phenomenon has been less frequently reported in Asia (Bataille 
et al, 2018; Scheele et al, 2019). Despite the widespread occurrence 
of Bd in Asia (Goka et al., 2009; Swei et al, 2011; Bai et al., 2012; 
Zhu et al. 2014, 2016), there remains a lack of observation of 
large amounts of individual deaths or amphibian population 
declines (Olson et al., 2013; Scheele et al., 2019). One potential 
explanation is that there may be certain unique community 
structures involved in amphibian skin symbiosis that contribute 
to low Bd load or decreased mortality. Therefore, exploring 
the characterization and investigation of the environmental 
correlations of the skin microbiota symbionts of native 
amphibians in Asian regions is crucial for understanding the 
natural resistance of these animals to Bd. 

There are still limited number reports on the skin 
microbiota symbionts of amphibians in Asia. In mainland 
China, previous studies have mainly focused on Rana dybowskii 
in northeastern China; this is an economically important frog 
species because its oviduct can be used in traditional Chinese 
medicine. These studies discussed the differences in skin 
symbiotic bacteria between farmed and wild individuals and 
seasonal changes in the gut and skin microbiotas of R. dybowskii 
(Bie et al., 2019; Tong et al., 2019). Zhao et al. (2019) cloned the 
16S rRNA gene to conduct bacterial community analysis of 
the skin of Odorrana grahami, a frog mainly distributed in 
the Yun-Gui Plateau. Studies on skin symbionts of frogs in 
other parts of China are still scarce, and most studies have 
focused on only one species. Previous research has shown that 
amphibian taxonomy (host species) is an important predictor 
of the diversity and community composition of skin bacteria 
(McKenzie et al, 2012; Kueneman et al, 2014; Bletz et al., 2017a), 
and geography (geographical regions), ontogeny (life histories) 
and ecomorphology (habitats) are also important predictors of 
the skin microbiome (Abarca et al., 2018; Longo et al., 2015; Bletz 
et al, 2017b). 

Here, we examined the skin bacteria of five sympatric 
amphibian species from Lishui, Zhejiang Province, a place 
with no alien amphibian invasions (Liu and Li, 2009; Liu et 
al, 2013) and a high diversity of amphibian species endemic 
to southeastern China. These five species have various habitat 
preferences, with three species, namely, Fejervarya multistriata, 
Hylarana guentheri and Microhyla fissipes, belonging to the 
aquatic type, generally inhabiting standing water, and one 
species, Polypedates megacephalus, belonging to the arboreal type, 
often living in low shrubs, grasses or agricultural crops and 
approaching standing waters during breeding seasons from 
April to June. The final species, Rana zhenhaiensis, is a terrestrial 
frog endemic to China, with adult animals inhabiting forest 
grasslands and using standing waters for reproduction only in 
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mating seasons from approximately January to April. Their 
diverse habitat use provides us with an ideal opportunity to 
investigate and compare the diversity of the skin bacterial 
communities of sympatric amphibian species, which will be 
helpful for future research on the relationship of the skin 
microbiome with Bd prevalence in southeastern China. 


2. Materials and Methods 


2.1. Sample collection Our study was conducted in April 
2019, and we randomly sampled individuals of each of five 
species (F. multistriata, H. guentheri, M. fissipes, P. megacephalus, 
R. zhenhaiensis) from Lishui, Zhejiang Province. The five species 
are mainly distributed in eastern and southern China and 
some countries in southeastern Asia. They showed low or no 
Bd prevalence and no evidence of Bd-induced mortality events 
in both historical specimens from museum and specimens 
from the wild (Bai et al., 2012; Zhu et al., 2014) (for sample 
information, see Table 1). Similar to F. multistriata, H. guentheri 
and M. fissipes, P. megacephalus was also sampled near water 
because the sampling time was its spawning season, while R. 
zhenhaiensis was sampled on a hillside, away from the water. 
We followed standard procedures and captured each individual 
by hand using new disposable polyethylene gloves and then 
poured ~100 mL of sterile purified water to remove dirt and 
transient bacteria from the surface of each frog (Belden et al., 
2015). After flushing, we swabbed each frog with a total of 30 
strokes, including five strokes on each side of the abdominal 
midline, the inner thighs and the foot webbings of each hind 
leg (Vredenburg et al., 2010). Swabs were stored individually in 
labelled sterile Eppendorf tubes in a car refrigerator before they 
were stored at -20 °C in the laboratory for further processing. 


2.2. Bacterial DNA extraction and sequencing Bacterial 
genomic DNA was extracted from swabs using the E.Z.N.A. 
Bacterial DNA Kit (catalogue No. D3350-02; Omega Bio-tek, 
Norcross, GA, USA) according to the manufacturer’s protocol 
with minor adjustments in the pretreatment of the samples. We 
first immersed the swabs in ImL of Ix TE buffer, vortexed and 
discarded the swabs, centrifuged the suspensions at 10000xg for 
3 min, and discarded the supernatant. Then, we added 100 pl 
of 1x TE buffer, vortexed the mixture to completely resuspend 
the pellet and followed the protocol provided by the kit. The 
extracted DNA was quantified by a NanoDrop (Thermo Fisher 
Scientific, Wilmington, DE, USA) and diluted to 5 ng/ul using 
RNase-free water. Then, the hypervariable V4 region of the 
16S rRNA gene of symbiotic bacteria was amplified using the 
barcoded primers 515F (5-GTGCCAGCMGCCGCGGTAA-3) 
and 806R (5-GGACTACHVGGGTWTCTAAT-3) (barcode- 
515F and barcode-806R) (Knutie et al., 2017). The final amplicon 
reaction contained 5x FastPfu Buffer (100 mM Tris-SO,, 200 
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mM KCI, 50 mM (NH,),SO,, 10 mM MgSO, 10% glycerol) 
(TransGen, Beijing, China), 2.5 mM dNTPs, 2.5 U/ul FastPfu 
DNA Polymerase (TransGen, Beijing, China), 20 ng of DNA, 
10 uM barcode primers and PCR grade water to 50 pl. The 
PCRs included a denaturing step at 95°C for 2 min; 35 cycles 
of 95 °C for 20 s, 51 °C for 20 s, and 72 °C for 30 s; and a final 
extension for 5 min at 72 °C. The concentration and purity of 
the amplification products were tested by a Qubit fluorometer 
and agarose gel electrophoresis, and then, the amplification 
products were purified and pooled in equimolar quantities and 
sequenced by Novogene Bioinformatics Technology Company 
(Beijing, China) on an Illumina HiSeq platform using a 250 bp 
paired-end strategy (Caporaso et al., 2012). Raw sequencing data 
are available in the NCBI Sequence Read Archive database 
with the BioProject accession number PRJNA613376. 


2.3. 16S sequences and statistical analyses After 
demultiplexing and trimming the barcodes and primers of 
each sample, the overlapping forward and reverse paired 
reads were assembled with Flash (Magoc and Salzberg, 2011). 
The assembled fastq files were processed using the open-source 
Quantitative Insights Into Microbial Ecology 2 (QIME 2 
(version 2018.11)) pipeline (https://qiime2.org; Hall and Beiko, 
2018; Bolyen et al. 2019). We used default parameters in QIME 
2 (q2-quality-filter plugin) to conduct quality control, and then, 
the sequences were denoised by Deblur in QIIME 2 (q2-deblur- 
denoise plugin) (Amir et al, 2017). Sequences shorter than 250 bp 
were discarded (Bletz et al, 2017c). Subsequently, the chimeras 
were filtered using UCHIME (Edgar et al., 2011). Then, the rest 
of sequences were assigned into the sequence variants using 
the Greengenes database (version 13_8) (McDonald et al., 2012; 
DeSantis et al., 2006) through the q2-feature-classifier plugin, 
and the taxonomic composition at the phylum, family and 
genus levels was generated based on operational taxonomic 
units (OTU) annotation. After removing the sequences derived 
from chloroplasts, mitochondria, archaea and eukaryotes (q2- 
taxa-filter plugin), we removed the OTUs that covered < 0.005% 
of the total reads (q2-feature-table-filter plugin) (Bokulich et al, 
2013). The number of reads per sample ranged from 21,420 to 
154,799, and all samples were rarefied to 21,420 reads to mitigate 
the effects of uneven sampling (Schloss et al., 2013). 

We used Faith’s phylogenetic diversity (PD) (Faith, 1992) 
index to construct the rarefaction curve by randomly 
subsampling the data at a series of sequence depths to test 
whether the depth was sufficient to capture the majority of taxa 
(Supplementary Figure S1). Before using the rarefied dataset to 
estimate alpha and beta diversity of the skin bacteria by QIME 
2, a phylogenetic tree was created for generating phylogenetic 
diversity measures such as unweighted and weighted UniFrac 
distances (Lozupone and Knight, 2005) or PD via the qiime 
alignment mafft, qiime alignment mask, qiime phylogeny 
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fasttree and giime phylogeny midpoint-root commands. Then, 
alpha and beta diversity metrics were assigned using the q2- 
diversity plugin. We estimated the observed OTUs, evenness, 
Shannon and PD of the skin bacteria to assess differences 
among samples. Kruskal-Wallis tests were used to compare 
alpha diversity among different species. Beta diversity 
was calculated with the Bray-Curtis, Jaccard, unweighted 
UniFrac (based on presence-absence and phylogenetic 
distance) and weighted UniFrac (based on abundance and 
phylogenetic distance) distance metrics and visualized with 
principal coordinates analysis (PCoA). We used permutational 
multivariate analyses of variance (PERMANOVA) (Anderson, 
2005) with 999 permutations to determine the role of the host 
species in shaping skin bacterial communities. We performed 
further analyses and data visualization in the R environment, 
version 3.6.1 (R Core Team, 2019). We used the R package 
“geplot2” (Wickham, 2016) to generate stacked bar charts at 
the phylum, family, and genus levels to visualize differences in 
the relative abundances of skin bacterial taxa. Using the same 
package, we visualized the differences in the alpha diversity of 
the five species. To examine differentially abundant bacterial 
taxa between the ectomorphs (habitat types) of samples, we 
used linear discriminant analysis (LDA) effect size (LEfSe) 
(Segata et al, 2011) on the Galaxy Web platform. We defined the 
dividing point as LDA 22.0 as recommended in the reference. 


3. Results 


We used 16S rRNA gene amplicon sequencing on the Illumina 
platform to examine skin bacterial communities from all 46 
swab samples of each individual. After sequence filtration and 
classification, we detected a total of 1020 OTUs belonging to 17 
phyla, with 9 phyla shared by all host species, 7 shared by 2-4 
species and one found on only M. fissipes. Among the phyla, 
Proteobacteria (56-87%), Firmicutes (7-27%), Actinobacteria (3- 
11%) and Bacteroidetes (2-5%) were the dominant skin symbiont 
bacteria of all five host species (Figure 1). At the family level, 
Enterobacteriaceae was the most abundant in M. fissipes, H. 
guentheri, F. multistriata, and R. zhenhaiensis (21-50%), whereas 
Exiguobacteraceae was the most prevalent in P. megacephalus 
(22%) (Figure 1). At the genus level, Escherichia was the most 
abundant in M. fissipes, H. guentheri, and F. multistriata; an 
unnamed genus in Enterobacteriaceae was the most abundant 
in R. zhenhaiensis (18-21%); and Exiguobacterium was the most 
prevalent in P. megacephalus (22%) (Figure 2). 

We observed significant differences in the alpha diversity of 
the skin bacteria among host species (Shannon, Kruskal-Wallis: 
H = 191905, P = 0.0007; PD, Kruskal-Wallis: H = 17.5626, P = 
0.0015; observed OTUs, Kruskal-Wallis: H = 19.2590, P = 0.0007; 
evenness, Kruskal-Wallis: H = 15.2537, P = 0.0042). Multiple 
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Table 1 Sampling locations for five amphibian species from Lishui, Zhejiang, China. 


Species Family N Coordinates 

Fejervarya multistriata Dicroglossidae 11 N28.463601° E119.936487° 

Hylarana guentheri Ranidae 11 N28.461126° E119.900026° 

Microhyla fissipes Microhylidae 9 N28.460355° E119.938849° 

Polypedates megacephalus Rhacophoridae 5 N28.439375° E119.913218° 

Rana zhenhaiensis Ranidae 10 N28.463867° E119.936191° 
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Figure 1 Relative abundance of skin bacterial taxa of five host species. A) Relative abundance of skin bacterial taxa at the phylum 
level across species. B) Relative abundance of skin bacterial taxa at the family level across species. Taxa with relative abundances < 1% 
were clustered together. 


385 


Xuejiao YANG et al. 


Characterization of Skin 
Microbiome of Frogs 


| No.4 


were highly abundant in F. multistriata, H. guentheri, M. fissipes 
and P. megacephalus samples from watersides areas; however, R. 
zhenhaiensis sampled from land contained more Proteobacteria 


than other species (Figure 5). 


2, Figure 3), and the results of the observed OTUs and evenness 
were similar to the Shannon index and PD results (Table S1). 


comparisons showed that the Shannon index and PD were 
significantly lower in R. zhenhaiensis than in other species (Table 


For beta diversity of the skin bacteria, PERMANOVA 


analyses showed that the host species indeed had a significant 
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on skin symbionts of five amphibian species in southeastern 
Bacteroidetes dominated all five host species, but the alpha and 


China. We found Proteobacteria, Firmicutes, Actinobacteria and 


among most species pairs, especially between R. zhenhaiensis and 


each of the other four species (Table 3, Table S2). 


beta diversity and the main taxa at the family and genus level 


Given the unique alpha and beta diversity of the skin 


of the skin bacteria were different among host species. Lower 
alpha diversity and richer Proteobacteria were found on the 
terrestrial species R. zhenhaiensis than on the four aquatic species, 
F. multistriata, H. guentheri, M. fissipes and P. megacephalus. 

This is among the first studies to seek the skin symbionts 
of multiple species and analyse them together in southeastern 


bacteria of R. zhenhaiensis, the only species sampled from land, 
we further explored the role of habitat type or sampling site 
type on the skin symbiont of the host using LEfSe to identify 
specific taxa that varied in abundance in the samples sourced 
from aquatic and terrestrial environments. Interestingly, 
Actinobacteria, Bacteroidetes, Cyanobacteria and Firmicutes 
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Figure 2 Heatmap of bacterial genera with relative abundances > 0.1% across five host species based on Bray-Curtis distance. 
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Table 2 Multiple comparisons of species pairs by Kruskal-Wallis tests for Shannon index and PD determination among host species. 


Boldface indicates significance < 0.05. 


| Vol. 11 


Index Group 1 Group 2 H p-value q-value 
Shannon F. multistriata (n=11) H. guentheri (n=11) 1.1739 0.2786 0.4643 
M. fissipes (n=9) 0.1169 0.7324 0.8633 
P. megacephalus (n=5) 0.0802 0.777 0.8633 
R. zhenhaiensis (n=10) 10.4926 0.0012 0.0048 
H. guentheri (n=11) M. fissipes (n=9) 1.3867 0.239 0.4643 
P. megacephalus (n=5) 0.0032 0.9548 0.9548 
R. zhenhaiensis (n=10) 14.4595 0.0001 0.0014 
M. fissipes (n=9) P. megacephalus (n=5) 0.36 0.5485 0.7836 
R. zhenhaiensis (n=10) 10.14 0.0015 0.0048 
P. megacephalus (n=5) R. zhenhaiensis (n=10) 6 0.0143 0.0358 
PD F. multistriata (n=11) H. guentheri (n=11) 0.0528 0.8182 0.9548 
M. fissipes (n=9) 0.6364 0.425 0.7084 
P. megacephalus (n=5) 0.0032 0.9548 0.9548 
R. zhenhaiensis (n=10) 10.4926 0.0012 0.006 
H. guentheri (n=11) M. fissipes (n=9) 0.7633 0.3823 0.7084 
P. megacephalus (n=5) 0.0032 0.9548 0.9548 
R. zhenhaiensis (n=10) 12.3967 0.0004 0.0043 
M. fissipes (n=9) P. megacephalus (n=5) 0.36 0.5485 0.7836 
R. zhenhaiensis (n=10) 8.1667 0.0043 0.0121 
P. megacephalus (n=5) R. zhenhaiensis (n=10) 7.935 0.0048 0.0121 
= Shannon = PD 
254 t25 
204 +20 
= 154 -15 
= — = 
E | | J 
Zw 10. | L10 
54 E] = e | L5 
a = | 


M. fissipes 


Species 


Figure 3 Shannon index and PD of skin bacteria among the five host species. 


H. guentheri P. megacephalus F multistriata R. zhenhaiensis 


China, where diverse amphibians are distributed. A detailed 
understanding of amphibian skin symbionts in different 
regions of the globe could aid the exploration of measures to 
mitigate chytridiomycosis (Abarca et al., 2018). Southeastern 
China has been proven to be a suitable habitat for Bd, as 
inferred by models (Liu et al., 2013); however, large-scale 
population declines and extinctions have not been detected 
by extensive field surveys (Bai et al, 2012; Zhu et al, 2014). One 


of the possible reasons for this is that Bd originated in Asia 


and has a long coevolutionary history with amphibians in 
Asia (O'Hanlon et al., 2018). Another important factor is the 
beneficial effect of the skin microbiome on the amphibian host 
through microbial species diversity, community composition 
and assemblage of anti-Bd bacterial taxa, which can produce 
anti-fungal compounds (Harris et al., 2006; Loudon et al., 
2014; Jani et al., 2017; Bell et al., 2018). Considering the close 
relationship between the skin symbionts and the host, it is of 
great significance to explore the skin symbionts of amphibians 
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Figure 4 Principal coordinate analysis (PCoA) plots of the beta diversity of the skin bacterial community. A) Bray-Curtis matrices 
and B) weighted UniFrac matrices of five host species. Each point represents the skin bacterial community of an individual sample. 


in Asia to understand host health and disease resistance in this 
region. 

It has been reported in various studies that host ecology 
or habitat and/or host taxonomy are crucial drivers of skin 
microbial composition and diversity in amphibians (McKenzie 
et al., 2012; Bletz et al., 2017b). We found clear differences in 
the skin bacteria of the five frog host species that we assessed 
and pronounced differences in ecomorphology between 
riparian and terrestrial species in terms of both of the skin 
microbial diversity and composition. With regard to the 
Shannon index and PD, the main difference was reflected in 
ecomorphology; however, contrary to previous studies that 
showed that terrestrial amphibians have greater bacterial 
diversity than arboreal or aquatic species (Abarca et al., 2018; 
Rebollar et al., 2019), we found that the bacterial diversity 


of R. zhenhaiensis was significantly lower than that of other 
riparian frogs. Nevertheless, the Shannon index of the skin 
bacteria of Craugastor fitzingeri, a terrestrial toad in Panama, 
exhibited a median value among the riparian, arboreal 
and terrestrial amphibians (Rebollar et al., 2016). One of the 
possible explanations for the results is that the environmental 
conditions, such as temperature, humidity and altitude, varied 
in the sampling areas. Alternatively, while our study and 
Rebollar et al. (2016) contained only one terrestrial species, 
the studies by Abarca et al. (2018) and Rebollar et al. (2019) 
lacked aquatic species. Therefore, studies on additional host 
species inhabiting different environments around the world 
are needed to test whether the difference in alpha diversity 
among ecomorphologies is ubiquitous or limited to specific 
species. Furthermore, unlike the results for the gut microbiota, 
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Table 3 Multiple comparisons by PERMANOVA of (-diversity among species based on Bray-Curtis and weighted UniFrac distance matrixes. 
Boldface indicates significance < 0.05. 


index Group 1 Group 2 Sample size pseudo-F p-value q-value 
Bray-Curtis F. multistriata H. guentheri 22 1.5416 0.042 0.07 
M. fissipes 20 0.4619 0.978 0.978 
P. megacephalus 16 P3552, 0.15 0.1875 
R. zhenhaiensis 21 3.2958 0.001 0.0033 
H. guentheri M. fissipes 20 12595 0.141 0.1875 
P. megacephalus 16 2.8489 0.001 0.0033 
R. zhenhaiensis 21 5.5809 0.001 0.0033 
M. fissipes P. megacephalus 14 1.0217 0.409 0.4544 
R. zhenhaiensis 1 2.5287 0.017 0.0425 
P. megacephalus R. zhenhaiensis 15 1.7858 0.04 0.07 
weighted UniFrac F. multistriata H. guentheri 22, 1.5872 0.105 0.15 
M. fissipes 20 0.3882 0.909 0.909 
P. megacephalus 16 1.1978 0.273 0.3413 
R. zhenhaiensis 21 5.3213 0.001 0.005 
H. guentheri M. fissipes 20 1.6926 0.034 0.0567 
P. megacephalus 16 2.8789 0.006 0.012 
R. zhenhaiensis 21 11.4543 0.001 0.005 
M. fissipes P. megacephalus 14 0.7413 0.572 0.6356 
R. zhenhaiensis 19 4.3736 0.005 0.012 
P. megacephalus R. zhenhaiensis 15 5.3313 0.002 0.0067 
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Figure 5 Cladograms showing that the bacterial taxa differ between species sampled from aquatic and terrestrial environments based 
on LEfSe. 


we found no evidence for strong phylosymbiosis in these phylogeny (Brooks et al, 2016). This is consistent with the result 
skin microbial communities of the five frog hosts, that is, we of Bletz et al., 2017b. In addition to the difference in coevolution 
did not find skin microbial communities parallel with host between gut and skin with the bacterial symbionts elucidated 
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by Bletz, it is plausible that skin bacteria varied constantly via 
environmental transmission (Loudon et al., 2014), and probable 
physical contact between species living in the same or adjacent 
habitats caused this effect. 

In the present study, we found 17 bacterial phyla on the 
amphibian skin, among which Proteobacteria, Firmicutes, 
Actinobacteria and Bacteroidetes were the four dominant 
phyla, which is consistent with observations in previous studies 
conducted in both temperate and tropical systems (eg, Walke 
et al, 2014; Belden et al, 2015; Bletz et al, 2017b; Jiménez et al, 
2019), suggesting that the skin of amphibians may be a selective 
niche for the colonization of peculiar bacterial taxa (Rebollar 
et al., 2016; Bletz et al, 2017b). However, the relative abundances 
of these phyla seemed to vary to some extent across different 
areas; Proteobacteria was almost the most abundant phylum 
on amphibians around the world, and the other main phyla 
may be shifted (Belden et al., 2015; Bletz et al, 2017b; Abarca et 
al., 2018). In this study, Proteobacteria was more abundant on 
the terrestrial species R. zhenhaiensis than on the other species, 
which is similar to the results of other studies (Belden et al., 
2015; Abarca et al., 2018). Proteobacteria was also the most 
abundant of summer and wild samples of skin microbiotas on 
the terrestrial species R. dybowskii in northeastern China (Bie 
et al, 2019; Tong et al., 2019). In southeastern China, Firmicutes 
was the second most abundant phylum detected on the skin 
of amphibians, which was different from the results obtained 
in some tropical countries, such as Panama, Costa Rica and 
Madagascar (Belden et al., 2015; Bletz et al, 2017b; Jiménez et al, 
2019). Apparently, the frogs of temperate regions may have 
higher Firmicutes abundance than those of tropical zones 
(Abarca et al., 2018). 

At the family and genus levels, Enterobacteriaceae 
and Exiguobacteraceae and the genera Escherichia and 
Exiguobacterium belonging to these two families were identified 
as the most abundant taxa on the skin of our focal species. 
This result is interesting in light of largely similar studies that 
showed that the most dominant taxon was Pseudomonadaceae 
or Pseudomonas in different areas (Becker et al., 2015; Rebollar et 
al., 2016, 2018; Prado-Irwin et al., 2017; Catenazzi et al., 2018; Bie et 
al, 2019). However, Enterobacteriaceae was the most abundant 
family on Mantella aurantiaca in Madagascar (Passos et al., 2018). 
Interestingly, Enterobacteriaceae and Pseudomonadaceae often 
exhibit high relative abundances together on the skin of many 
amphibians distributed in many countries, such as Madagascar, 
Germany, USA, and Australia (Bletz et al, 2017b, c, d; Kearns 
et al, 2017; Bell et al, 2018). Additionally, a high proportion 
of inhibitory isolates belonging to Enterobacteriaceae and 
Pseudomonadaceae was found in the database of culturable 
anti-Bd bacteria (Woodhams et al., 2015). With regard to the 
other abundant genus, Exiguobacterium, in our study, also 
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showed high relative abundance in a population of common 
midwife toads (Alytes obstetricans) (Bates et al., 2018). To date, a 
total of 17 species have been found to belong to Exiguobacterium, 
and the isolates of this genus have been found to exhibit 
bioremediation and enzyme production and degradation 
of toxic substances. Exiguobacterium is a versatile bacterial 
genus distributed in a variety of environments (Kasana and 
Pandey, 2018). The high abundance of Enterobacteriaceae and 
Exiguobacteraceae in this study may indicate the putatively 
unique anti-Bd function among the five sympatric frog species 
in southeastern China. Of course, further study is needed to 
verify this hypothesis through bacterial isolation and culture 
techniques. 

To our knowledge, this is the first study to assess differences 
in the skin microbiota of cooccurring species in southeastern 
China by high-throughput sequencing. It is important to 
describe and understand indigenous skin bacterial communities 
to develop conservation strategies that can be applied locally. 
Further sampling across a wide area of China is needed to 
promote the understanding of amphibian skin microbial 
ecology and to obtain an improved understanding of the 
particular mechanism of interaction between Bd and host 
amphibians in China. 
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Appendix 


— Fejervarya multistriata — Hylarana guentherit —*- Microhyla fissipes 
— Polypedates megacephalus ~*~ Rana zhenhaiensis 
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Figure S1 Rarefaction curves based on the phylogenetic diversity measure of skin bacterial samples from five host species. 
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Figure S2 Phylogenetic tree constructed by Phylogenetic Generalized Least Squares (PGLS) processed by R based on the amphibian 
tree constructed by Peloso et al., 2015. We have not found the species Microhyla fissipes in Peloso et al., 2015, therefore we used the 
closely related species Microhyla heymonsi to instead of Microhyla fissipes (Matsui et al., 2011). 


Table S1 Multiple comparisons of species pairs by Kruskal-Wallis tests for observed OTUs and evenness among host species. Boldface 
indicates significance < 0.05. 


Index Group 1 Group 2 H p-value q-value 


M. fissipes (n=9) 0.417 0.5184 0.632 


R. zhenhaiensis (n=10) 9.3953 0.0022 0.0054 


P. megacephalus (n=5) 3.4941 0.0616 0.1232 


M. fissipes (n=9) P. megacephalus (n=5) 1.2844 0.2571 0.4285 
p 


P. megacephalus (n=5) R. zhenhaiensis (n=10) 9.3918 0.0022 0.0054 


M. fissipes (n=9) 0.2439 0.6214 0.8647 


R. zhenhaiensis (n=10) 9.6 0.0019 0.0097 


P. megacephalus (n=5) 0.0032 0.9548 0.9548 


M. fissipes (n=9) P. megacephalus (n=5) 0.04 0.8415 0.935 


P. megacephalus (n=5) R. zhenhaiensis (n=10) 4.86 0.0275 0.0687 


Table S2 Multiple comparisons by PERMANOVA of {-diversity among species based on Jaccard and unweighted UniFrac distance matrixes. 
Boldface indicates significance < 0.05. 
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P. megacephalus R. zhenhaiensis 2.3715 0.017 0.045 


